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Abstract 


A numerical simulation of a convective turbulence 
event is investigated and compared with observational 
data. The specific case was encountered during one of 
NASA ’s flight tests and was characterized by severe 
turbulence. The event was associated with 
overshooting convective turrets that contained low to 
moderate radar reflectivity. Model comparisons with 
observations are quite favorable. Turbulence hazard 
metrics are proposed and applied to the numerical 
data set. Issues such as adequate grid size are 
examined. 

Introduction 

A major portion of the accidents from aircraft 
turbulence encounters are within close proximity of 
atmospheric convection (e.g., thunderstorms). 1 The 
National Aeronautics and Space Administration 
(NASA), through its Aviation Safety Program, is 
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testing technologies that will reduce the risk of injuries 
from these types of encounters. Primary focus of the 
turbulence element within this program is the 
characterization of turbulence and its environment, and 
the development and testing of hazard-estimation 
algorithms for both radar and in situ detection. The 
ultimate goal is to operationally test onboard sensors 
that will provide ample warning prior to encounters 
with hazardous turbulence. In support of turbulence 
characterization, numerical modeling of atmospheric 
convection is being conducted using a large eddy 
simulation model. A special need for the modeling is 
to provide realistic data sets for developing and testing 
turbulence detection sensors. However, the first step 
prior to this application must be verification that the 
numerical model can produce useful and realistic data 
sets. To meet this goal numerical simulations with 
NASA’s Terminal Area Simulation System (TASS) are 
being applied to actual cases with moderate and severe 
turbulence encounters. Validation of the numerical 
simulations relies on measurements from radars, 
satellites, and aircraft penetrations. 

One particular case targeted for numerical study is 
a turbulence event encountered during NASA’s fall- 
2000 flight tests. In this event, hazardous turbulence 
was encountered by NASA Langley’s B-757 on 14 
December 2000. Severe intensities of turbulence were 
measured as the aircraft penetrated updraft plumes 
near the top of a narrow line of thunderstorms. Data 
from onboard Doppler radar, in situ wind and 
temperature measurements, and recorded NEXRAD 
radar data (Tallahassee (TLH)) are available for 
comparison with the numerical simulation of this case. 
This paper will describe the conditions associated with 
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the turbulence encounter and analyze results from the 
simulation of this event. 

Description of Event 

During two test days in December 2000, regions of 
convective turbulence were purposefully encountered 
by NASA’s B-757. Areas with moderate or greater 
radar reflectivity, i.e. RRF > 35 dBZ , were avoided as 
routinely done by commercial air carriers. Turbulence 
measurements from the in situ system were quantified 
in terms of RMS normal loads (o ng ), where 0.20 g <<J ng 
<0.30 g is considered moderate and a ng > 0.30 g is 
severe. During two flights, 14 moderate and 4 severe 
turbulence encounters occurred in the vicinity of deep 
convection. Further details of these encounters can be 
found in a companion paper presented at this 
conference. 2 

The turbulence event selected for study was the 
strongest event encountered during the flight tests. 
This event, 191-6 (referred as 191 .3 in reference [2]) , 
had a peak turbulence intensity of c ng = 0.44 g, which 
is comparable to that in several incidents involving 
commercial aircraft (see Fig. 15 in Hamilton and 
Proctor 2 ). The in situ measurements of vertical wind 
and RMS normal loads during this encounter are 
shown in Fig. 1 . The turbulence was characterized by 
sharp oscillations in vertical velocity over a distance of 
about 5 km. 

Event 191-6 was associated with a narrow line of 
convective cells extending east-northeast (ENE) from 
the Gulf of Mexico through the Florida Panhandle (Fig. 
2). Storm tops reached 1 1 .8 (39,000 ft) mean sea 

level (MSL) with cell movement from the west- 


In Situ of Event 191-6 



Figure 1 . Measured vertical wind and RMS normal 
load acceleration vs distance along the flight path 
during event 191-6. 


southwest (WSW) at about 20 m/s. Turbulence was 
encountered as the NASA B-757 flew through weak 
radar reflectivity regions near the top of the convective 
line. At the time of penetration, the aircraft was at an 
altitude of 10.3 km above ground level (AGL), and was 
at a location just north of the Florida-Georgia state line 
(latitude, longitude: -84.17 W, 30.76 N). The aircraft 
heading was ENE (nearly parallel to the line) with an 
air speed of about 235 m/s. Further details of the flight 
test and meteorology surrounding this event are in 
Hamilton and Proctor 2 . 

The approach taken in this paper is to perform a 
numerical simulation of the convective line and 
examine subsequent turbulence fields near the altitude 
of the 191-6 encounter. 

The Model and Initial Conditions 

The numerical simulation is carried out with the 
Terminal Area Simulation System 3 ’ 4 (TASS), which is a 
large eddy simulation model developed for simulating 
convective clouds and atmospheric turbulence. 

Model Equations 

The TASS model consists of 12 prognostic 
equations: three equations for momentum, one 

equation each for pressure deviation and potential 
temperature, six coupled equations for continuity of 
water substance (water vapor, cloud droplet water, 
cloud ice crystals, rain, snow and hail) and a prognostic 
equation for a massless tracer. The model also 
contains numerous microphysics models for computing 
cloud and precipitation physical interactions. Further 
details of the model formulation can be found in 
references [ 3, 5, and 6 ]. 



near the time of event 191-6. 
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The TASS equation set (ignoring Coriolis terms) in 
standard tensor notation is as follows: 


Momentum: 
dm + H_ dp 
& Po dxi 
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dx 
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For precipitating variables (such as rain and snow), an 
additional vertical flux term is added to account for fall 
out. 


Buoyancy Term: 
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Thermodynamic Equation (Potential Temperature): 
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with the Potential Temperature defined as: 

Poo p 

p 


g = t 
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In the above equations, w, is the tensor component o: 
velocity, t is time, p is deviation from atmospheric 
pressure P, T is atmospheric temperature, p is the ail 
density, C p and C v are the specific heats of air at constam 
pressure and volume, g is the earth's gravitationa 
acceleration, R d is the gas constant for dry air, P 00 is £ 
constant equivalent to 1000 millibars (10 5 pascals) oi 
pressure, Q v is the mixing ratio for water vapor, 0 r is sun: 
of the mixing ratios for liquid and ice substances, L is the 
latent heat, and S is a water substance source term. 
Environmental state variables, e.g., p m Q vm P 0 , and 6 0 , 
are defined from the initial input sounding and are 
functions of height only. 

Conservation of Scalar Variables (e.g., water vapor, 
cloud droplet water, etc.): 


A modified Smagorinsky first-order closure is used 
for the subgrid eddy viscosity as: 


Km Is 
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The subgrid eddy viscosity for momentum, K Mt is 
modified by the Richardson numbers, for stratification, 
Ri Si and for flow rotation, Ri n with = 3 and a 2 = 1 .5. 
The subgrid eddy viscosity for heat and water substance 
is determined as K H = 3 K M . The subgrid length scale, l s 
is determined from the grid volume. 

Numerical Approximations 

Time-derivative approximations for momentum and 
pressure are time-split explicit for computational 
efficiency. The prognostic equations are approximated 
using 4 th -order energy-conserving space differencing 
and 2nd order time differencing. Only light numerical 
filtering is applied using a 6th-order filter. Potential 
temperature and water substances equations are 
approximated with third-order accurate time and space 
differences with upstream-biased quadratic interpolation. 7 
The horizontal derivatives in TASS are approximated on 
an Arakawa-C grid. 8 The numerical formulation for 
TASS is stable for long-term integrations and is 
essentially free from numerical diffusion. 9 


Model Domain and Boundary Conditions 

The domain is rotated 66° clockwise such that the 
x-coordinate is orthogonal to the convective line and 
they-coordinate is parallel to the line (see Fig. 3). The 
physical size of the domain is 25 x 25 x 14 km. The 
grid size is uniform at 1 00 m over most of the domain, 
except below an altitude of 2000 m where grid 
stretching shrinks the vertical size to 50 m. The 
domain is resolved by 148 vertical levels having 25 1 x 
251 grid points. 
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Figure 3. Orientation of TASS domain relative to 
convective line. Observed composite radar reflectivity 
field and flight path depicted near the time of event 
191-6. 



Figure 4. Skew-T chart with atmospheric sounding 
for time and location near event 191-6. Generated 
from forecast with Mesoscale Atmospheric Simulation 
System. Sounding data provide by Mike Kaplan, 
NCSU, under NASA contract NAS1-99074. 

Periodic boundary conditions are assumed at the 
northeastern and southwestern boundaries (orthogonal 
to the convective line), while open nonreflective 
conditions are assumed for the northwestern and 
southeastern boundary. 



Figure 5. Cumulative distribution of normalized 
aircraft loads as a function of wave number. Aircraft 
calculation based on B-7 57-200 frequency domain 
model at an altitude of 6.1 km with an airspeed of 212 
m/s and a weight of 180, 000 lbs. Assumes von Karman 
turbulence spectrum with an outer scale of 200 m and 
o w =l. 

The ground boundary is impermeable with nonslip 
velocity specifications. The surface stress due to the 
ground is determined locally from the wind speed, 
surface roughness, and the local thermal stratification. 
Details of the surface formulation are in the appendix of 
Proctor and Han. 10 

Initial Conditions 

The simulation is initialized with a vertical 
distribution of temperature, dewpoint, and wind velocity, 
representative of the environment near the time and 
location of the turbulence event. Since observed profiles 
were unavailable at the time and location of event 191-6, 
a forecast sounding (Fig. 4) was obtained from a 
mesoscale prediction model. 11 

Subgrid Fusion 

In order to model scales of motion important for 
aircraft response, high resolution is needed. Results 
from a frequency-domain flight dynamics model, 12,13 
indicated that scales of motion as small as 50 m (wave 
number of 0. 1 26 rad/m) are needed in order to capture 
at least 97% of the cumulative aircraft load distribution 
(Fig. 5). Available computer capability and the size of 
the computations restricts the grid size to about 100 m. 
Although this resolution misses scales that are 
important for aircraft response, the model’s ability to 
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simulate the larger-scale features of a convective- 
turbulence event can be assessed. This evaluation is 
discussed in the first part of the next section. High- 
resolution turbulence fields are achieved by extracting 
a subdomain from the numerical simulation and 
merging it with subgrid turbulence fields 14 . Results 
from this procedure, including onboard radar and 
aircraft simulations, are discussed in the second part of 
the next section. 

Results 

Simulation of Convective Line 

Results from the TASS 100 m simulation of the 
convective line are presented below, and are compared 
with observations derived from ground based 
NEXRAD radar and B-757 flight data. Table 1 shows 
comparisons between model and observations of 
selected features. 


Table 1. Model Comparison 


Variable 

TASS 

Observed 

Orientation of Line 

WSW- 

ENE 

WSW- 

ENE 

Peak Storm Tops 

W.5 km 

11.8 Am 

Peak Radar Reflectivity 
at Ground 

53.5 dBZ 

55 dBZ 

Peak Radar Reflectivity 
atz=9 km 

38.9 dBZ 

40 dBZ 

Cell Motion (toward) 

ENEdX 
19 m/s 

ENEzX 
17 m/s 

Width of Convective 
Line near Ground Level 

6 km 

8 km 

(based on 20 dBZ) 



Peak Eddy Dissipation 
Rate (m 2/3 /s) 

0.86 

0.74 

Horizontal Scale of 
Turbulence Patch at 

5 km 

5-6 km 

Flight Level 




Comparisons between modeled radar reflectivity 
and that observed by the Tallahassee radar are depicted 
in Figs 6-9. The model values are without consideration 
of the radar beam width, pulse length, and beam tilt. 
Nevertheless, only small differences exist, and 
comparisons show similar orientation, scale, and 
magnitude. 

The radar reflectivity from the onboard radar just 
before encountering event 191-6 is shown in Fig. 10. A 


simulation of this radar using the TASS data set is 
shown for comparison in Fig. 1 1 . Both show similar 
scale and intensity, although details in the echo 
structure differ. 

A three-dimensional perspective of the simulated 
convective line appears quite realistic, exhibiting 
cumulus turrets, anvil outflow and overshooting tops 
(Fig. 12). The convective cells exhibit downwind tilt 
(toward the northeast) with most of the anvil outflow 
spreading in that direction. During the actual 
encounter, the NASA B-757 flew toward the northeast 
parallel to the line and entered the overshooting cloud 
areas near the storm tops. Severe turbulence was 
encountered as the aircraft skirted the northwestern 
flank of the convective line. 

In Fig. 13, a horizontal cross section of vertical 
velocity is show at the flight level of the NASA B-757. 
Also superimposed is the TASS radar reflectivity 
(truncated at 15 dBZ). The radar reflectivity is 
primarily due to the presence of snow, with 
temperatures at this altitude being colder than -40 °C. 
Localized regions of strong upward velocity (peak of 
17 m/s) are associated with rising convective turrets 
and are embedded within low-radar reflectivity regions. 
Strong downward motions are especially noted on the 
downwind side of the turrets. This feature also was 
observed during the flight test and is discussed in 
Hamilton and Proctor. 2 Radar reflectivity within the 
downdraft regions is weak in comparison with the 
updraft regions. Also note that strong gradients of 
vertical velocity may occur in regions where radar 
reflectivity is less than 20-15 dBZ. Similarly, the in 
situ vertical velocity derived from the B-757 (Fig. 1), 
also shows intense pulses of updraft and downdraft 
during this turbulence event. 

Figure 14 shows TASS eddy-dissipation rate 
(EDR) for the same region as depicted in Fig. 13. 
Although the strongest values are located within higher 
reflectivity regions (i.e., greater 25 dBZ), moderate to 
strong values may occur in regions of weak radar 
reflectivity. As indicated in Table 1 , the peak values of 
simulated EDR are similar to those obtained from in 
situ data of the actual encounter. 

An energy spectrum computed from the TASS 
velocity data (Fig. 1 5) indicates a continuous spectra of 
turbulence rather than an isolated gust. This also is 
confirmed from spectrum of flight data that was 
measured during the actual event. 2 The TASS spectra 
(Fig. 15) appears to have an inertial subrange with a 
-5/3 slope especially at larger wavelengths. At smaller 
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TASS Radar Reflectivity at 
T = 49 min and z= 9 km 



-15 -10 


10 


Figure 7. Same as Fig. 6, but at 9000 m altitude. 



wavelengths, however, the spectra shows a steeper 
slope than the theoretical -5/3 slope. This drop-off in 
energy at higher wavenumbers is often found in other 
LES studies, 15,16 and is theoretically expected since 
values at each grid cell represent volumetric averages 
rather than point values. 17 

Merging with Subgrid Turbulence 

Although the 100 m TASS simulation was able to 
simulate the larger-scale features of the turbulence 
event, it could not resolve the smaller-scales of motion 
important for aircraft response calculations. Figure 15 


shows that only wavelengths greater than 600 m (6 grid 
points) are adequately resolved, and according to Fig. 5 
only 40% of the cumulative aircraft load is captured at 
these frequencies. Since finer resolution is needed for 
proper aircraft response simulation and turbulence 
radar simulation, high-resolution subgrid turbulence 
fields were merged with a subdomain of the TASS 
simulation. This data set was generated by NCAR 
using the following procedure: 18 

1) A sub-volume of the domain was selected 
which encompassed the turbulence event. 
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2) The variables (at t=48 min) were interpolated 
to a 25 m grid, within a 12.8 km x 12.8 km 
horizontal and 3.2 km vertical subdomain. 

3) Following a technique devised by Frehlich et 
al 14 , subgrid wind fields using a von Karman 
algorithm were then merged with the TASS 
data. The von Karman subgrid parameters 
(variance and outer length scale) were 
determined from a best fit of the model 
generated structure functions (after 
interpolation) to the desired Kolmogorov 
behavior. 


As shown in Fig. 16, the merging of the subgrid 
extends the inertial subrange to scales less than 50 m. 
The data set is now sufficient to resolve most of the 
scales important for aircraft response. 

A horizontal cross section of the vertical velocity 
field from the merged data set is shown in Figs. 17. 
With the addition of the subgrid windfields, the peak 
updraft and downdraft speeds at flight altitude are 27 
m/s , and -17 m/s, respectively. Significant fluctuations 
of vertical velocity are confined within regions having 
some precipitation or cloud material, (cf. Fig. 1 8). 


Subgrid fields are only added to the velocity fields. 
Other fields, such as radar reflectivity are simply 
interpolated to the higher-resolution subdomain. 



-W- 

Pm-DBZ 


Figure 10. Radar reflectivity (dBZ) from onboard 
turbulence radar. Observed just prior to encounter 
with event 191-6. Range rings every 4 bn. Image 
provide by Les Britt, RTI, under NASA contract NAS1- 
99074. 


Hazard Analysis 

Using algorithms developed by Bowles, a hazard 
metric for aircraft turbulence can be applied to the data 
set. For a particular aircraft, the RMS normal load can 
be estimated from g w using look-up tables; 13 i.e., 

g( x >y) = F{c w , altitude, aircraft type, weight, 
airspeed} 

The g w fields can be computed for any horizontal plane 
in the merged data set, by using a moving average as: 


°w( x >y) = 


1 2 2 __ 

, T J J {w(x',y')-w(x,y)} 1 2 dx'dy' 

L * L y L x L v 

x — - v — y_ 

2 y 2 



Figure 11. Radar reflectivity (dBZ) from simulation 
of onboard radar using TASS data set. Simulation 
assumes same altitude and heading as in Fig. 10. 
Range rings every 2 km. Image provide by Les Britt, 
RTI, under NASA contract NAS 1-99 07 4. 


where the averaging interval along the jc and y 
coordinates is L x , L y , respectively. The average 
vertical wind, w , is computed from the vertical wind, 

w, as: 

l x l 

X+— yH — — 

1 r 2 2 

Mx,y) = ~— J J w{x\y')dx'dy' 

2 7 2 

The value for the averaging interval, L x =L y = 
1000 m, is chosen to correspond to a 4-5 ^ averaging 
period for a commercial jet at cruise speeds. Hence, 
the second moment of the w-field is computed 
assuming a 1 x 1 km moving box. 

The g w field (Fig. 19) computed from w in Fig. 
17 exhibits a peak value of 8.2 m/s. 
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Figure 12. TASS generated convective-cloud line for event 191-6 as viewed from southeast. 


The RMS normal load (o ng ) is computed from the 
G w fields, assuming aircraft parameters for NASA’s B- 
757 (Fig. 20). Since calculations are independent of 
aircraft heading, evaluation of the turbulence field is 
relatively simple. Regions with G ng > 0.3 g represent 
severe turbulence, while 0.20 g <G ng <0.30 g represent 
moderate turbulence. The peak value from the data set 
is 0.37. The relationship of the hazardous regions with 
the radar reflectivity can be seen by comparing Figs. 1 8 
and 20. Note that moderate intensity of turbulence is 
sometimes found in regions of very weak radar 
reflectivity. 

Flight Simulation 

A one-dimensional profile of vertical velocity 
(Fig. 21) was extracted along the arrow shown in Fig. 
18. The profile has a peak updraft of 24 m/s and is 
surrounded by strong gradients in vertical velocity. 
The data in Fig. 21 was used as input into a time- 
domain flight dynamics model. The output from the 
dynamic model is shown in Fig. 22. The dynamic 
model gives a peak G ng = 0.363 which compares 
favorably with the peak value of 0.37 shown in Fig. 20. 
Very close agreement is obtained even though 
independent methods are used in calculating a ng . This 
supports the credibility of the technique described 
above for computing G ng from numerical-model data. 
A comparison of Fig. 22 with the observed profile of 
Gng (Fig.l) exhibits interesting similarities, although 
weaker than the measured G ng = 0.44. 

Also of interest is the relation between the peak 
aircraft loads and its RMS value (see Fig. 22). The 


peak aircraft loads are nearly instantaneous; and 
individually, would be difficult to detect by look-ahead 
sensors. In contrast, the g^ predictions of hazard 
encompass the regions of strong peak loads, and are 
broad enough in horizontal scale to be easily detected. 

Radar Simulation 

Radar simulations using the merged data set are 
currently underway. A preliminary result of radar 
spectral width is shown in Fig. 23. Spectral width 
actually observed during event 191-6 are shown for 
onboard-turbulence radar and TLH ground-based radar 
in Figs. 13b and 14b of reference [2]. The peak value 
in the radar simulation (7 m/s) compares favorably with 
that observed by TLH radar (7 m/s) and onboard radar 
(8-9 m/s). 


Table . 2. RMS Normal Load Comparison 


Source 

Peak a ne (g’s) 

In situ 

0.44 

Onboard Turbulence Radar + 

0.37 

Ground-Based Doppler Radar + 

0.33 

Flight Dynamics Simulation 

0.36 

Model Diagnostic from g w field 

0.37 

Radar Simulation with Model 

0.33 

Data 



+ Computed from radar spectrum width (see reference 13 ) 


RMS Normal Load Comparison 

Table 2 shows a comparison between peak RMS 
normal loads measured near the B-757 flight path with 
those simulated from the merged data set. Values for 
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Figure 13. TASS radar reflectivity (flood) and vertical 
velocity (2 m/s contours, negative values dashed) 
along a horizontal cross section at flight level. 


Spectra: TASS Simulation of R-191-6, A=100 m 



X (meters) 

Figure 15. Energy spectra from TASS, as computed 
from velocity field within 25x25 km horizontal plain at 
flight elevation (t=49min). 



Figure 14. TASS eddy dissipation rate to the 1/3 
power (m 2/3 s f at time and location corresponding to 


Fig 13. 

G n g may be obtained from the radar spectrum width by 
using look-up tables developed by Bowles. 13 All 
sources indicate a severe turbulence event whether 
from observed data or simulation. 


Summary and Conclusions 


A numerical simulation of a convective turbulence 
event is investigated and compared with observational 
data. The results have been validated with data from 
ground based NEXRAD radar, onboard flight radar, 
and in situ measurements. The numerical results show 


Fir 6-1 91 Spectra at Elevation 10,300 m 
merged subgrid 



10 * 1 0‘ 2 10 -’ 

Wavenumber 


Figure 16. Energy spectrum from vertical velocity 
field shown in Fig. 1 7. 

severe turbulence associated with buoyant plumes that 
penetrate the upper-level thunderstorm outflow. The 
peak radar reflectivity in these plumes is 36 dBZ. The 
simulations produce updraft plumes of similar scale to 
those encountered during the test flight. The simulated 
radar reflectivity compares well with that obtained 
from the aircraft’s onboard radar. 

Resolved scales of motion as small as 50 m are 
needed in order to accurately diagnose aircraft normal 
load accelerations. Given this requirement, realistic 
turbulence fields may be created by merging subgrid- 
scales of turbulence to a convective-cloud simulation. 
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A hazard algorithm for use with model data sets is 
demonstrated. The algorithm diagnoses the RMS 
normal loads from second moments of the vertical 
velocity field and is independent of aircraft motion. 



Figure 17. Vertical velocity field from merged data 
set. Horizontal cross section at flight altitude (z= 10.3 
km AGL). [Line shows a hypothetical flight path.] 
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1 Km Averaged Second Moment of Vertical Velocity 


TASS Simulation at 48 minutes with merged subgrid 




Figure 19. Diagnosed o w for time and domain shown 
in Fig. 17. 



Figure 18. TASS radar reflectivity field (dBZ) 
interpolated to domain shown in Fig. 1 7) 



Figure 20. Turbulence hazard field (<j ng ) computed 
from G w in Fig. 19. 
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Figure 21. Vertical velocity field extracted from 
merged data set along flight path shown in Fig. 17. 



Figure 22. Normal load acceleration and RMS 
normal load acceleration ( (J ng ) vs distance along the 
flight path as computed from the vertical velocity 
shown in Fig. 21. Calculation assumes 5 sec moving 
average (~ 1 km) and aircraft parameters equivalent to 
that in 191-6. Local calculations based on 2 degree of 
freedom B-757 time-domain flight dynamics model. 



R-Mm (m) = 1075. Alt. (ft) 33B00. 


SPEC TRA L WIDTH DiSPLA Y 

Figure 23 Spectrum width as simulated with B-757 
onboard radar. Range rings every 2 km. Image 
provide by Carol Kelly, RTI, under NASA contract 
NAS1-99074. 
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